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Abstract 

Numerical  rootfinding  problems  are  quite  common  in  stochastic  modeling. 
However,  many  solutions  stop  at  the  presentation  of  a  probability  generat¬ 
ing  function  for  the  state  probabilities.  But  with  increasing  easy  access  to 
computing  power,  many  problems  whose  answers  were  typically  left  in  in¬ 
complete  form  or  for  which  there  has  been  a  search  for  alternative  solution 
methods  are  currently  being  reexamined.  The  class  of  Markov  chains  whose 
transition  matrices  have  quasi-triangular  layouts  (i.e.,  those  having  sub-  or 
super-triangular  sets  of  zeros)  is  a  good  case  in  point.  They  have  an  espe¬ 
cially  nice  structure  which  leads  to  a  rather  concise  representation  for  the 
generating  functions.  But  the  complete  solution  then  requires  the  finding  of 
roots.  Fortunately,  these  problems  can  be  shown  to  have  special  properties 
that  make  accurate  roothnding  quite  feasible,  and  we  thus  supply  an  efficient 
numerical  procedure  for  solution.  1 
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1  INTRODUCTION 


As  numerous  authors  have  noted  (for  example,  Abolnikov  and  Dukhovny, 
1987,  Bailey,  1954,  and  Powell,  1985),  many  denumerable  discrete-time  Markov 
chains  (with  particular  applications  in  inventory,  dam  and  queueing  model¬ 
ing)  have  one  of  two  special  transition-matrix  structures.  These  forms  have 
been  typically  called  quasi-triangular  because  of  the  presence  of  sub-  or  super- 
triangular  sets  of  zeros: 
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S,  =  1  U  brl. 

r>  < 1 

The  structure  of  these  matrices  loads  1»>  some  paiticulaily  r«>miso  mp- 
resentations  for  the  probability  generating  functions  (PGFs)  of  the  Markov 
chain  equilibrium  state  process.  When  t  he  stationary  equation  for  such  a 


ivlarkov  chain  is  exercised,  the  PGF  of  the  steady-state  probabilities  has  an 
algebraic  function  in  its  denominator  whose  roots  are  critical  in  the  final  so¬ 
lution.  (We  henceforth  refer  to  this  denominator  function  equated  to  zero  as 
the  system’s  characteristic  equation). 

For  the  transition  matrix  A,  the  characteristic  equation  (CE)  turns  out 
to  be  (at  least  for  complex  z  with  absolute  value  <  1) 

OO 

a':'  =  a(2)>  (*) 

1=0 

where  K  is  as  defined  in  the  matrix  representation  A  (corresponding,  for  ex¬ 
ample,  to  a  constant  batch-input  module  in  bulk  queues).  Under  the  assump 
tion  that  a(z)  possesses  all  its  derivatives  at  z  =  1  (i.e.,  that  the  distribution 
{a,}  has  all  moments),  a(z )  may  be  set  equal  to  the  Laplace-Stieltjes  trans¬ 
form  of  a  distribution  function  [call  that  A(t),  and  set  its  mean  to  1/p]  of  a 
nonnegative  random  variable  evaluated  at  A(1  —  z),  where  A  is  an  arbitrary 
positive  constant  for  the  time  being.  Thus  we  may  also  write  that 

zK  =  A'\\{\  -  z)\.  (2) 


The  CE  associated  with  the  matrix  B  may  be  written  as 

OO 

ZK  =  Tbtz'  -  0(z),  (3) 


zK  =  B' j/i(]  -  z)\.  (4) 


In  these  representations,  the  constant.  K  is  as  given  in  B  (corresponding, 
for  example,  to  a  constant  batch-service  module  in  bulk  queues),  0{z )  is 
defined  as  the  PGF  of  the  probabilities  {fi, }.  and  B+  is  the  Laplace-Sl  ielt  jes 
transform  of  a  distribution  function  [(all  t  hat  B(t),  and  set  its  meant"  1  A1  nf 
a  nonnegative  random  variable  evaluated  at  //( 1  z),  where  //  is  an  arbitrary 

positive  constant. 


Recognize  that  Equations  (1)  and  (3)  are  generalizations  of  the  well- 
known  fundamental  equation  of  branching  processes,  typically  written  as 
z  —  f(z ),  where  /  would  be  the  PGF  for  the  number  of  offspring  emanating 
from  one  parent.  Gross  and  Harris  (1985).  for  example,  provide  the  details 
of  the  root  problem  for  this  model.  In  actuality,  it  is  the  B  problem  which 
is  in  fact  the  more  direct  relative  of  the  branching  process,  so  it  is  this  one 
on  which  we  comment  in  detail  first.. 

The  CE  of  the  matrix  B  may  be  rewritten  in  the  standard  way  by  using 
z  =  r  exp(i6),  and  we  find  that 

rK exp(iOK)  =  £?*[/r(  1  -  r  exp(id))}exp(2irni)  (5) 

for  n  —  1,2 . K .  This  equation  clearly  has  a  root  at  unity,  and  by  Rouche's 

theorem,  we  can  show  that  there  are  K  others  inside  the  unit  circle  j  z  |—  1 
when  the  chain  is  ergodic.  The  condition  for  ergodicity  is  that 

oo 

(3'{  1)  -  Y  nb"  >  K 

n-  ( 1 


or 

-  * 

These  can  be  shown  to  be  equivalent  to  the  requirement  that  K A / p  <  1. 
For  the  A-matrix  problem,  recall  that  the  characteristic  equation  is 

=  A*[A(i  -  z)J,  (6) 

where  A*  is  a  Laplace-Stieltjes  transform.  Ergodocity  obtains  here  when 

Q,(  1 )  —  Y  7?an  <  K 
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This  is  equivalent  to  requiring  that  A  A’//  •  1. 


It  is  easily  shown  by  Rouche’s  theorem  that  (6)  has  K  roots  inside  and 
on  the  unit  circle,  including  the  root  z—  1.  Abolnikov  and  Dukhovny  (1987) 
have  noted  that  all  the  roots  on  the  unit  circle  are,  in  fact,  simple. 

In  the  prior  work  of  Cnaudhry,  Harris  and  Marchal  (1989),  we  have  seen 
under  an  assumption  of  infinite  divisibility  that  Equation  (1)  from  matrix  A 
mav  be  rewritten  as 


2 


K 


M:)]* 


It  then  follows  that 


z^a1(z)e„iK  (n  =  1.2, A'),  (7) 

where  e„  K  is  an  n La  (out  of  K)  complex  rc>ot  of  unity. 

When  the  original  a(2)  is  infinitely  divisible,  the  function  afiz)  is  itself 
a  legitimate  PGF  with  at(l)  =  A//V/a  <'  1.  The  equation  2  =  <*1(2)  has  only 
the  root  2  =  1;  the  remaining  K  —  1  roots  inside  and  on  the  unit  circle  follow 
distinctly  from  the  other  K  -  1  equations  using  Rouche’s  theorem. 

The  analyticity  of  the  Kth  root  of  0(2)  and  subsequent  use  of  Rouche’s 
theorem  for  distinctness  also  follows  when  a(z)  is  nonzero  though  not  neces¬ 
sarily  infinitely  divisible.  (See  Chaudhry,  Harris  and  Marchal,  1989.) 

Unfortunately,  it  is  not  true  that  all  0(2)  associated  with  such  models 
have  no  zeros  inside  the  unit  circle.  However,  we  can  feel  comfortable  know¬ 
ing  that  a  good  number  of  the  problems  encountered  in  practice  will  have 
infinitely  divisible  distributions  since  many  such  probability  functions  are 
built  up  by  mixtures  and  convolutions  from  the  infinitely  divisible  exponen¬ 
tial  and  Erlang.  For  all  other  distributions,  one  should  always  first  try  to 
determine  whether  a(z )  is  ever  zero,  for,  if  not,  the  Kth-root  approach  of  the 
infinitely  divisible  case  will  work. 

The  vanishing  of  0(2)  is  equivalent  to  requiring  that 

TOO 

A*[X(l-z)]=  c~Xii-»dA(t)  =  0. 

J  n 


By  changing  to  the  poiar  form  2  =  ?■( cos 0  4  isind)  and  then  separating  ma! 
and  imaginary  parts,  this  is  identical  t<>  asking  that,  simultaneously. 
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Jn 


-  A(1  -r  co*  0)t 


cos|(  Ar  sin  ft)t\ dA(i)  —  H 


and 


Mi-rco ’<»tsln\(Xr  sm9)t}a(t)dt  -  0. 


Chaudhry,  Harris  and  Marchal  (1989)  have  established  firm  sufficient  con¬ 
ditions  for  distinctness  of  the  roots  for  ( 1 )  and  (2).  But  neither  of  their  con¬ 
ditions  is  necessary.  Examples  are  presented  in  their  paper  for  which  a(z) 
has  zeros,  but  where  the  characteristic  equations  still  have  distinct  roots.  In 
fact,  they  never  found  an  example  without  distinct  roots.  But  a  proof  that 
this  is  indeed  always  true  is  elusive,  and  the  search  for  a  necessary  condition 
will  be  the  subject  of  future  research. 

All  of  these  arguments  hold  in  a  parallel  way  for  the  13-matrix  problem, 
as  well,  as  long  as  /3'(1)  >  A',  with  its  A  roots  strictly  inside  the  unit  circle. 
The  B-matrix  version  of  the  problem  would  then  have  (3\(z)  in  place  of  ai(z) 
in  Equation  (7). 


2  ALGORITHM 

The  Chaudhry,  Harris  and  Marchal  (1989)  experiences  for  both  the  A-  and 
B-matrix  problems  yielded  a  nearly  circular  pattern  of  roots  on  the  inside 
of  the  unit  circle.  This  is,  indeed,  perfectly  predictable  given  the  above 
results.  This  comes  directly  from  writing  the  right-hand  side  of  the  Kth 
root  of  the  characteristic  equation  as  the  product  of  a  root  of  unity  and  the 
Kth  root  of  the  probability  generating  function.  The  support  of  the  random 
l!?  rorrptpAr d’^g  to  the  n<rf  r> ,  is  increasingly  focused  towrard  0,  so  that 
the  net  effect  of  the  multiplication  is  to  move  the  root  slightly  away  from 
the  unit  circle’s  boundary  toward  the  renter,  largely  preserving  a  circle-like 
appearance.  In  examples  offered  in  the  following,  we  display  a  number  of 
diagrams  of  such  patterns.  Some  possible  guidance  for  the  speeding  up  of 
the  algorithm  cornes  to  mind  from  this  formulation.  It  is  to  use  the  nth  root 
(out  of  K,  n  —  1 ,2,...,K)  of  unity  as  the  starting  point  for  each  attempted 
solution  of  (7).  Thus  we  would  set  the  initial  value  of  the  radius  at  1  and  0 
at  ImKjn.  We  can  even  do  a  bit  better  since  we  know  that  modulus  of 
each  root  must,  be  less  than  the  positive  real  root;  therefore  we  us*'  an  initial 
modulus  for  any  particular  problem  equal  to  the  smallest  one  obtained  so  far 
as  n  moves  to  K. 


Compare  this  to  the  approach  used  by  Chaudhrv  and  colleagues  i  jo> 
example,  see  Rriere  and  Chaudhrv  (  WWi  to  find  the  roots  for  the  (7  1 

problem,  as  in  Equations  (31  and  H  r  Their  algorithm  works  on  the  logarithm 
of  Equation  (4),  using  the  real  and  imaginary  parts  separately  First,  the 
real  equation  is  solved  for  r,  holding  ft  fixed  Then  the  imaginary  equation 
is  solved  for  n,  given  the  derived  r  and  fixed  ft.  It  the  result  is  integral,  we 
are  finished,  otherwise,  the  fractional  part  of  r>  is  used  and  compared  to  its 
counterpart  related  to  a  larger  ft.  If  the  former  is  smaller,  then  9  is  incieased 
further.  Such  comparisons  are  repeated  until  there  is  a  decrease  between  the 
fractional  parts  of  the  last  two  derived  n  values.  This  indicates  that  a  root 
is  located  at  an  angle  between  the  last  two  9  values.  Details  of  the  algorithm 
are  given  in  Briere  and  Chaudhrv  t  10o7  r 

We  have  built  up  our  algorithm,  then,  from  the  intuitively  appealing 
near-circular  pattern  found  for  the  roots  The  roots  are  ordered  accoiding  to 
increasing  angle  6.  with  a  starting  value  foi  the  search  for  the  (n  H)st  root 
using  the  modulus  of  the  nth  root  and  the  angle  of  the  (n  +  l).st  root  of  unity. 
fn-i,K-  We  chose  to  explore  the  use  of  a  fixed-point  algorithm,  knowing 
that  a  sufficiently  well-chosen  starting  point  would  lead  to  convergence  of 
successive  substitutions  of  z  into  Equation  (7).  This  is  particularly  easy  to 
formulate  in  light  of  the  linearity  of  the  left-hand  side  of  (7). 

Experimentation  with  this  fixed-point  technique  has  led  to  success,  and 
detailed  examples  are  presented  in  the  following  sections.  To  start  the  overall 
algorithm,  it  is  necessary  to  provide  one  of  the  roots  of  the  problem  and  the 
obvious  init ial  choice  is  to  find  the  real  root  on  (0,  1 ),  since  we  know  the  exact 
angle  of  this  particular  root.  Any  method  could  be  used  to  solve  for  this  real 
root,  actually,  but  we  chose  to  use  the  same  fixed-point  approach,  simply  for 
uniformity  in  the  algorithm.  The  derivation  of  a  safe  starting  point  for  the 
location  of  the  real  root  proved  to  be  an  interesting  problem  in  its  own  right. 
Due  to  the  knowledge  that  the  support  for  O)  is  toward  0,  an  initial  guess 
greater  than  the  value  of  the  root  would  be  natural. 

At.  first,  an  initial  value  equal  to  for  the  A-matrix  problem  (or 

equivalently.  K!8'(\)  for  the  B-matrix  model;  seemed  logical,  sir."'  this  is 
the  (single)  real  root  ( --  p)  for  the  basic  hi  A/.M  queueing  problem  While 
this  starting  value  did  work  in  ail  e;;r.'-s.  Dm  accuracy  of  the  real  r< » >1  h >’iml 
was  not  ideal  because  small  but  both*  rs-mie  complex  i omp- ■nrr.ts  crop'  into 
the  solution.  Nevertheless,  all  successive  root.s  wore  found  correctly  ami 
their  values  matched  those  obtained  via  other  root  finding  packages  (r.g..  the 


^t’ACK  implement  at  inn  of  the  Chaudhrv  approach) 

The  progress  of  convergence  of  the  iteration  from  staiting  values  outside 
of  the  near-circular  pattern  of  roots  could  suffer  from  discontinuities  in  the 
derivatives  of  the  function  «i(r)  brought  about  by  the  zeros  (if  any)  of  the 
function  a.  But  it  has  been  observed  that  the  roots  of  at  z )  lie  always  outside 
of  the  locus  of  the  roots  of  the  characteristic  equation  (7).  It  has  also  been 
observed  that  the  root  locus  was  skewed  asymmetrically  toward  the  positive 
axis  in  the  complex  plane,  thus  making  the  roots  with  angles  closest  to  rr 
those  of  shortest  modulus.  Thus,  if  one  wished  to  provide  starting  values 
within  the  basin  of  attraction  of  the  desired  roots,  but  avoid  possible  dis¬ 
continuities  in  the  derivative,  one  could  seek  starting  values  at  the  angle  of 
2  "  2KitJ  1  K/21  radians,  finding  the  roots  in  the  top  half  of  the  complex 

plane  by  seeking  the  roots  in  the  order  2 rr I\  v  with  n  ■—  I\  2,  h  2  -  1 . <*. 

thus  ending  with  the  real  root  in  (0.  1  ). 

This  variation  on  the  technique  proved  particularly  successful,  and  a  start¬ 
ing  value  of  (r.B)  —  (0,  —  9K:2)  for  the  first  root  was  used.  The  effect  of  this 
variation  was  to  present  a  value  of  slight ly  less  modulus  than  the  desired 
root,  and  at  an  angle  of  the  corresponding  root  of  unity,  as  the  starting  value 
for  the  fixed  point  iteration  on  each  root..  In  general,  the  numerical  impact,  of 
the  variation  was  to  improve  the  effi c  iencv  of  convergence  for  the  fixed  point 
by  one  to  two  iterations  for  each  complex  root.  The  real  root  on  (0,1]  was 
found  much  more  quickly  by  “coming  from  the  inside"  of  the  root  locus,  and 
the  complex  component  of  the  root  was  exactly  zero. 

3  Proof  of  Convergence 

For  purposes  of  this  proof,  we  assume  that  a{z)  is  infinitely  divisible  and 
thus  that  any  one  of  its  roots  is  a  proper  PGF. 

Theorem  1  Define  the  nth  characteristic  equation  for  the  G/E^.-l  and 
E k  / (7 /  1  problems  as 
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hi  :v  i 


or 


at  z„ ' 


i  h 


/' 


°  I  <  ~r,  I'r,  h- 


where  enK  is  the  nth  of  K  roots  n f  unity,  ami  the  roots  zKl.z^2  j,.  ,Z\.z.,. 
K2  =  [A/2J.  are  ordered  starting  noth  the  root  of  greatest  angle  m  the  upper 
complex  half  plane,  progressing  toward  iht  positive  real  axis.  From  an  initial 
point  ;<0)  =  (r<u8v)  with  modulus  r,,  «•  :c„,|h  1  A  n  4  i  <  A2.  the  method 

of  successive  substitutions, 

*i,M)  ■=  ,  'M 

applied  to  the  nth  characteristic  equation  for  the  (!  A’/.,  !  and  A’/.,-  G  1  prob¬ 

lems  vjiII  converge  to  a  unique  root. 

Proof:  Realize  first  that  the  characteristic  equation  is  separable  into  the 
factors  a'  and  (say'.  Note  that  ai(z„)  is  a  function  of 

zifl  at  the  i th  iteration  of  successive  substitution,  but  that  en  h-  is  unchanged 

by  the  iteration.  We  are  assured  that  all  the  zn,  n  -  1 . A  are  distinct 

simple  roots  in  or  on  the  unit  circle,  as  previously  proven  in  Ohaudhry.  Harris 
and  Marchal’s  Theorems  1  and  4.  1  It  thus  suffices  to  show  then  that  the 
fixed-point  approach  is  guaranteed  to  converge  to  a  root  in  jrj  ■'  !,  since  the 
distinctness  of  t.he  roots  will  assure  us  that  there  cannot  bp  convergence  to 
the  same  answer  for  more  than  one  en  !<.  Also  note  that  it  has  been  shown 
empirically  that  |zn+i[  <  \zn\  for  fl  <  n  f  A’2 

To  complete  the  proof,  then,  the  critical  need  is  to  show  that  the  proba¬ 
bility  generating  function  ct-i(z)  is  less  than  1  -  S  (f>  an  infinitesimal  positive 
constant)  for  any  2  chosen  within  or  on  the  closed  and  bounded  contour 
\zl.  <  1  -  6.  This  is  indeed  true  because 

x  ai.«z‘j  -  x  i~l  x  °«.*(i  - 

1=0  !  t-o  t-o 

oc 

<  (1  -  6)  X,  Q1.1  =  1  — 

t-n 

since  the  {alt}  define  a  legitimate  probability  distribution  under  infinite 
divisibility  whenever  the  traffic  intensity  <.  1.  OO 

’For  certain  0(2).  this  mav  be  due  to  the  fact  that  the  function  0(2)  is  nonzero  over  Hu 
complex  unit  disk  (at  least),  and  for  rrianv  problems  because  o ( e )  is  infiirii-lv  dr.isibh 
In  an v  ot  her  case,  i ,e.,  if  n(z )  itself  has  roots  111  t  hp  complex  plane,  then  the  derivatives 
f(z)  ar»  not  continuous  everywhere.  This  alone  does  not  preclude  convergence  to  a  fixed 
point,  but  may  result  in  erratic  performance  of  the  iteration. 


4  RESULTS 


For  the  exercise  of  our  algorithm  (called  Rl'IX),  we  chose  to  apply  the  tech¬ 
nique  to  a  variety  of  the  systems  analyzed  in  Chaudhry,  Harris  and  Marcha! 
(198u)  -  CHM  in  the  fo'lowing.  Th**  problems  covered  here  u.re  listed  by 
the  corresponding  table  in  that  r ferei^ed  paper:  Table  1  (APK^/M/l),  Ta¬ 
ble  2  {PH/Ek/  1),  Table  3  {GH/Er/ 1),  Table  5  {PH/EK/ 1),  and  Table  7 
(Rn/Exf  1),  all  solved  as  G/Ek/  1  or  B-matrix  formulations;  and  Table  8 
{Ek/GE/V)  solved  as  an  Er/G'I  or  A-matrix  model.  We  show  our  results 
in  tables  numbered  identically.  To  the  problems  of  CHM,  several  El,  Ej< ,  1 
and  H2 !  Ek  i  1  problems  unique  to  this  paper  were  added  to  illustrate  the 
exercise  of  our  algorithm  further.  1  hose  additional  problems  appear  in  this 
paper  as  Tables  4  and  6,  respectively  For  a  basis  of  comparison,  the  algo¬ 
rithm  of  Chaudhry  and  colleagues  was  exercised  on  the  same  systems  as  our 
fixed-point  algorithm. 

The  stopping  rule  for  the  iteration  on  each  root  was  that  successive  root 
values  be  less  than  10  ~15  in  the  square  of  Euclidean  distance  between  them 
in  the  complex  plane.  This  resulted  in  approximately  a  10~7  error  in  the  real 
and  complex  comnonents  of  the  results 

In  addition,  an  implementation  of  a  bivariate  Newton’s  method  due  to 
Kahaner,  Moler  and  Nash  (1989)  was  run  over  all  cases.  The  times  shown 
are  those  for  computing  and  displaying  the  roots  in  the  upper  half-plane. 

To  provide  a  clearer  picture  of  the  environments  under  which  our  results 
were  obtained,  we  indicate  in  Table  0  the  nature  of  the  hardware  used  for 
each  algorithm.  The  contents  of  Table  tl  apply  to  all  problem  types,  except 
for  Tables  IV  and  VI,  where  the  Chaudhry  algorithm  was  run  on  a  16  mHz 
386  without  a  coprocessor.  Thus  the  rootfinding  times  in  Tables  IV  and  VI 
are  without  any  bias  from  differing  computational  platforms. 


Tabic  0: 

Computational  Platforms 

RFIX  CHACO  KM&N 
machine  16  mHz  12  mHz  16  mHz 
on  which  386  286  386 

algorithm  no  with  no 

was  run  coprocessor  coprocessor  coprocessor 

Table  I: 

Selected  Scenarios  from  Table  1,  CHM,  M/Ex/l 
(A  =  .l) 


A 

K 

intensity 

KXJ  fi 

RFIX 

(in  sec) 

CHAUD 
(in  sec) 

KM&N 
(in  sec) 

20.0000 

10 

.05 

1 

39 

5 

50.0000 

25 

2 

43 

12 

200.0000 

100 

7 

51 

51 

1000.0000 

500 

30 

953 

251 

1.0526 

10 

.95 

2 

4 

7 

2.6316 

25 

3 

4 

15* 

10.5263 

100 

7 

10 

54* 

52.6316 

500 

30 

281 

257* 

Found  wrong  positive,  real  mol  loss  than  1.  This  happened  despite  an 
al  guess  ~  mod(z2)  <  zx . 


Table  II 

Selected  Scenarios  from  Table  2,  CHM,  PH j E^j\ 


(A: 

:  4.21) 

K 

intensity 

KXffi 

RF1X 
(tn  sec) 

UIAUD 

( m  sec) 

KM&N 

(in  sec) 

.460 

10 

.91522 

3 

21 

6* 

4.210 

10 

.10000 

3 

21 

5 

.31 T 

15 

.90214 

4 

30 

9' 

2.806 

15 

.10000 

2 

28 

7 

.155 

30 

.90214 

5 

49 

16* 

1.403 

30 

.10000 

3 

49 

14 

f  Found  real  root  =  1.0  instead  of  zj  <  z0  =  1. 


Table  III: 

Selected  Scenarios  from  Table  3,  CHM,  GH/E^/l 
(A  =  .85714) 


M 

A 

intensity 

KX/n 

RF1X 

( in  sec) 

CHAUD 
( in  sec) 

KM&N 
(in  sec) 

9.00 

10 

.95238 

5 

21 

5# 

85.00 

10 

.10084 

1 

21 

4 

14.01 

15 

.91771 

4 

52* 

6* 

128.00 

15 

.10045 

2 

29 

6 

29.00 

30 

.88670 

4 

go- 

12 

257.00 

30 

.10006 

3 

48 

14 

*  Algorithm  missed  some  complex  roots,  which  were  captured  in  a  second 
pass. 

#  Positive,  real  root  not  found  (error  message  =  ITERATION  NOT 
MAKING  GOOD  PROGRESS). 


Table  IV: 

Selected  Scenarios  for  Model  Flzj Fy  /  1 
(A  .1) 


K 

intensity 

A'A/p 

RFIX 

(in  see) 

CHAUD# 

(in  sec) 

KM&N 

(m  sec ) 

2.8571 

20 

.7 

3 

5 

20’ 

6.4286 

45 

5 

8 

51* 

14.2857 

100 

10 

13 

117* 

20.0000 

20 

.1 

2 

25 

23* 

45.0000 

45 

4 

28 

62* 

100.0000 

100 

9 

33 

90* 

2.2222 

20 

.9 

3 

5 

22* 

5.0000 

45 

6 

7 

48* 

11.1111 

100 

10 

12 

123* 

+  KM&N  became  totally  lost,  finding  several  roots  two  or  three  times, 
despite  initial  guesses  at  unique  values.  The  KM&N  subroutines  found  the 
negative,  real  root  (a  second  time)  when  given  a  positive,  real  initial  value,  for 
example.  Generally,  KM&N  did  all  right  in  the  third  quarter-plane  (negative¬ 
negative),  becoming  confused  at  or  before  8  —  90°  for  K=20  or  8  —  45°  for 
K  =  100.  Then  previously  found  roots  would  be  duplicated,  with  the  negative 
real  root  (if  any)  recurring  most  often. 

#  The  Chaudhry  software  was  run  "ii  a  16  mHz  386  without  a  c<>pp«es- 


Table  V: 

Selected  Scenarios  from  Table  5,  CHM,  PH /Ex/ 1 
(A  =  0.2) 


A 

K 

intensity 

RFIX 

CHA1JD* 

KM&N# 

KX/fi 

(in  sec) 

(in  sec ) 

(in  sec) 

20.0 

10 

.100 

2 

25 

2.2 

10 

.909 

5 

26 

30.0 

15 

.100 

2 

35 

3.3 

15 

.909 

5 

38 

60.0 

30 

TOO 

3 

63 

6.6 

30 

.900 

6 

66 

*  QPACK  solution  was  for  slightly  different  parameter  values  using  the 
same  model.  The  resultant  root  values  were  not  identical  to  those  obtained 
using  RFIX  and  KM&N,  but  we  feel  that  the  above  QPACK  run  times  are 
representative. 

*  KM&N  could  not  perform  the  rootfinding  operation  for  this  problem 
type.  Incorrect  values  were  found  for  all  points. 


Table  VI: 

Selected  Scenarios  for  Model  H2/Ek/1 
(Aj  =  X2  —  •  1 ) 


K 

intensity 

RFIX 

CHAUD* 

KM&N 

KX/fi 

( in  sec) 

(in  sec) 

(in  sec) 

1.0 

15 

.100 

3 

27 

10 

1.0 

30 

.100 

5 

29 

19 

1.0 

95 

.100 

14 

37 

60 

1.0 

500 

.100 

57 

660 

0.2 

15 

.500 

4 

7 

10 

0.2 

30 

.500 

6 

8 

21 

0.2 

95 

.500 

14 

17 

63 

0.2 

500 

.500 

57 

333 

.11 

15 

.909 

5 

5 

11* 

.11 

30 

.909 

7 

7 

21* 

.11 

95 

.909 

15 

74* 

64* 

.11 

500 

.909 

58 

365# 

*  For  the  real  root  on  the  uppersheet  at  360°,  there  was  agreement  be¬ 
tween  KM&N  and  RFIX  only  to  four  or  so  digits. 

#  Missed  real  roots,  so  routine  was  rerun  with  smaller  step  size. 

|  The  Chaudhry  software  was  run  on  a  16  mHz  386  without  a  coprocessor. 


I ' 


Table  VII: 

Selected  Scenarios  tor  Mode!  Rnj &k/  1 
(A  -  2) 


A 

A' 

intensity 
K  X/n 

RFIX 
(in  stcj 

CHAUD 

(in  sit) 

KM&N 

(in  sec) 

20.00 

10 

.10000 

1 

19 

5 

2.20 

10 

.90909 

2 

20 

6 

20.00 

15 

.10000 

2 

26 

7 

2.13 

15 

.93750 

2 

27 

8 

20.00 

30 

.10000 

3 

46 

14 

2.06 

30 

.96774 

3 

45 

15 

2.20 

100 

.90909 

8 

131 

46 

Table  VIII: 

Selected  Scenarios  from  Table  8,  CHM,  Ex  j GjB3/1 

(/'  -  2) 


\ 

K 

intensity 

X/h'fx 

RT  IX 

(in  sec ) 

CHAUD 
(in  see ) 

KM&N 
m  sec) 

2 

10 

.1 

2 

23 

8 

18 

10 

.9 

2 

22 

n 

1 

3 

15 

.1 

2 

31 

9 

27 

15 

.9 

2 

31 

9 

6 

30 

.1 

3 

51 

20 

36 

30 

.9 

3 

50 

21 

5  OBSERVATIONS 


The  performance  of  the  fixed-point  algorithm  can  be  seen  to  be  quite  effec¬ 
tive  and  stable  accross  a  wide  range  of  intensities  and  root  quantities.  A 
unique  feature  of  the  fixed-point  algorithm  is  that  the  number  of  successive 
substitution  iterations  required  per  root  actually  goes  down  as  the  number 
of  roots  grows  (i.e.,  as  K  grows).  This  is  understandable  since  the  differences 
between  the  moduli  decrease  as  K  grows,  and  the  starting  point  for  each  root 
is  closer  to  the  final  point. 

The  performance  in  this  problem  of  the  other  algorithms,  namely,  CHAUD/ 
QPACK  and  KM&N,  is  indicative  of  a  fundamental  difference  between  equa¬ 
tions  typically  requiring  rootfinding  in  applied  probability.  In  the  QPACK 
and  KM&N  implementations,  the  equation  zh  -  a(z)  is  analyzed  for  zeros 
by  tuAh.g  whether  0  =  u(z)  -  zK .  This  equation  represents  the  usual  ge¬ 
ometric  conceptualization  in  the  complex  plane  with  K  zeros  inside  and  on 
the  unit  disk.  Contrast  this  with  the  equation  used  by  RF1X, 


=  V~r 


—  '!K 


for  n  —  0,1,...,  A'  —  1.  Recall,  as  stated  in  Sec’ ion  3,  that  once  the  root 
of  unity  en,K  is  specified,  zn  =  a: (zn)f-n,K  is  one  of  the  K  branches  of  the 
multiple- valued  function  a{zn)x^K .  (For  example,  see  Churchill,  Brown  and 
Verhey,  1974,  p.  88.)  Each  of  the  K  branches  yields  an  equation  of  the  form 
zn  ~  Q]  (zn)  c n.K,  each  of  which  is  a  mapping  of  the  domain 
{r  >  0,  --7T  <  9  <  7r }  onto  the  domain 

{  yT  —  p  >  0,  (2n  —  1  )n/K  <  4>  ■  (2v  4  1  )tt / A’},  n  =  0, ...,  A’  -  1 . 

Since  each  of  these  latter  domains  ran  be  transformed  by  an  angular  ro¬ 
tation  to  {p  >  0,—ir/K  <  <t>  <  4  7t  /\  }.  the  RFIX  algorithm  amounts  to 
solving  K  completely  separate  problems  with  a  single  root.  This  is  a  geo¬ 
metric  interpretation  of  the  proof  in  Section  3  that  the  algorithm  will  always 


converge  to  the  intended  root,  whereas  the  other  algorithms  are,  in  fact,  solv¬ 
ing  for  the  K  roots  of  a  single  problem,  thus  allowing  for  the  possibility  of 
“skipped  roots”  and  convergence  to  undesired  roots. 

To  close  out  our  computational  work,  we  offer  pictures  of  the  precise 
locations  of  the  roots  for  three  of  the  problems  included  in  Table  VIII,  namely, 
Figure  1  for  K  —  10,p  =  .1;  Figure  2  for  K  =  15,p  —  .9;  and  Figure  3  for 
K  =  30,  p  =  .9. 


Figure  2: 


Plot  of  Roots  Inside  and  On  the  Unit  disk 
For  Eis/G E}/ 1  Model. 
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Figure  3: 


Plot  of  Roots  Inside  and  On  the  Unit  disk 
For  f.Kj/C£3/l  Model 


6  CONCLUDING  REMARKS 


The  efficacy  of  root  finding  for  this  problem  should  help  remove  the  impres¬ 
sion  that  rootfinding  in  stochastic  analysis  was  frought  with  hidden  obstacles. 
Results  have  indeed  been  uniformly  favorable. 

In  the  process  of  this  work,  it  was  observed  over  many  problems  that  all 
three  algorithms  spent  significant  time  searching  for  the  real  root(s)  of  the 
problem.  In  fact,  as  noted  in  Table  8,  RFIX  exhibited  convergence  prob¬ 
lems  for  queueing  systems  of  the  type.  In  general,  the  order  of 

convergence  of  fixed-point  iterations  rs  know  to  be  linear  (for  example,  see 
Nonweiler  and  Horwood,  1984),  while  it  is  known  quadratic  for  Newton’s 
method  and  1.84  (nearly  quadratic)  for  Muller’s  method. 

However,  our  computational  timings  suggest  a  much  better  generalk  record 
of  performance  for  our  fixed-point  approach.  As  pointed  out  in  Nonweiler 
and  Horwood  (1984),  the  order  o(  convergence  only  relates  to  errors  and 
values.  But  each  method  (namely,  fixed  point,  Newton’s  and  Muller’s)  has 
an  associated  computational  burden  per  iteration,  and  the  simplicity  of  the 
fixed-point  iteration  of  RFIX  is  shown  in  the  illustrative  code  in  Appendix  1 . 
By  contrast,  the  number  of  operations  in  the  high-level  language  implemen¬ 
tations  of  ^M&N’s  software  (Newton’s  method)  and  of  Chaudhry’s  QPACK 
(Muller's  method)  likely  far  exceed  those  of  RFIX  per  iteration;  indeed,  the 
sheer  sixe  of  the  KM&N  source  code  (comments  excluded)  can  give  some  clue 
as  to  the  fact  that  execution-time  instruction  count  is  several  times  greater. 

What  our  results  show  empirically  is: 

(a)  the  number  of  execution-time  instructions  per  iteration  is  far  more 
dominant  than  the  number  of  iterations,  producing  the  computational  com¬ 
plexity  relationship  Fixed  Point  <  N ewton' s  M ethod  <  Muller1  s  Method', 
and 

(b)  the  ratio  of  improvement  in  speed  of  convergence  (that  is,  timej  be¬ 

tween  methods  is  not  a  constant  from  problem  to  problem,  particularly  with 
changes  in  K .  That  is,  if,  for  a  fixed  K  of  any  problem  type,  we  look  at  a 
row  of  the  appropriate  table  (namely,  I  -  VIII)  and  calculate  Cj  and  c2  such 
that  t/tri.x  -  C)  *  In™, tor,1,  and  inns  -  c2  *  then  we  will  find  that 

over  all  I\  in  a  given  problem  (that  is.  table)  that  cj  and  c2  would  not  be 
constant. 

Thus,  over  all  I\  ,  we  might  say  that  t urix  —  cR/fi )*t^,wtor,<a  and  tnrix 
c2(l\)  *  where  cR/V)  and  r2(  f\  I  are  nonlinear  functions.  Thus  the 


computational  stability  (that  is,  consistency  of  speed  of  convergence)  of  our 
fixed-point  method  is  substantially  “more  than  linear"  relative  to  the  other 
methods  in  this  context. 

The  fact  that  this  feature  was  observed  over  a  wide  variety  of  character¬ 
istic  equations  implies  that  the  fixed-point  methdod  may  be  more  favorable 
than  other  approaches  for  larger  classes  of  rootfinding  problems  in  probability 
and  statistics. 
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Appendix  1 

Illustrative  FORTRAN  Code  for  Successive  Substitution  Algorithm 


program  algor 
complex  z,e,zl 
real  el,e2 
r0=.99 
pi=3. 1415923 
do  2  n  =  0,10 
t0=n*pi/10.0 
a=rO*cos(tQ) 
b=rO+sin(tO) 
z=cmplx(a,b) 
el  =  cos(tO) 
e2=sin(t0) 
e=cmplx(el,e2) 
do  1  i=l,100 

zl=e*(H(20.0/7.0)*(l-z))**(-.5) 

write(l,*),zl 

al=real(zl) 

a=real(z) 

bl  =  aimag(zl ) 

b  =  aimag(z) 

if(z.eq.zl)then 

goto  3 

else 

z^zl 

endif 

1  continue 

3  r0=sqrt(al**2+bl**2) 

write(  1  ’ 

2  continue 
stop 

end 
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